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Abstract. We investigate the stability properties and the dynamics of Bose-Einstein 
condensates with axial symmetry, especially with dipolar long-range interaction, using 
both simulations on grids and variational calculations. We present an extended 
variational ansatz which is applicable for axial symmetry and show that this ansatz 
can reproduce the lowest eigenfrequencies of the Bogoliubov spectrum, and also the 
corresponding eigenfunctions. Our variational ansatz is capable of describing the 
roton instability of pancake-shaped dipolar condensates for arbitrary angular momenta. 
After investigating the linear regime we apply the ansatz to determine the dynamics 
and show how the angular collapse is correctly described within the variational 
framework. 



PACS numbers: 03.75.Kk, 67.85.De, 05.45.-a 
1. Introduction 

As an alternative to the direct numerical solution of the Gross-Pitaevskii equation (GPE) 
and the Bogoliubov-de Gennes equations (BDGE), which describe the ground states 
and excitations of Bose-Einstein condensates (BECs) [l], Ran et al (2|[3] proposed a 
variational approach with coupled Gaussians to calculate stationary solutions and the 
stability of a condensate using the methods of nonlinear dynamics. Specifically, they 
used the time-dependent variational principle (TDVP) to map the GPE to a nonlinear 
dynamical system for the variational parameters whose fixed points correspond to 
ground states of the GPE. The stability can be analysed by means of the linearised 
equations of motion, i. e. in terms of the eigenvalues of the Jacobian. However, it 
remained unclear whether there is a direct relationship between the eigenvalues of the 
BDGE and those of the Jacobian. In a recent work [4] we analysed this problem 
for spherically symmetric systems and showed that a variational ansatz with coupled 
Gaussians has limitations in that it can only describe excitations with an angular 
momentum I < 2. We present an extended variational ansatz and show that there is 
indeed a good agreement between the lowest eigenvalues of the BDGE and the Jacobian 
for arbitrary angular momenta, with or without additional long-range interaction (LRI). 

A variational ansatz with a single Gaussian |5] cannot reproduce the biconcave 
structure of dipolar condensates revealed in certain parameter regions (6ll7j. Using 
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the variational ansatz of coupled Gaussians, however, it was possible to obtain non- 
Gaussians shapes (2]j3][8| and to reproduce the wave function and the energy of structured 
ground states. Additionally, the stability was analysed which revealed that the collapse 
breaks the axial symmetry within the variational approach. However, the ansatz used 
in Ref. js] is only capable of describing modes with angular momentum m = and 
m = 2 (see discussion in |4|). The angular collapse with m = 3 symmetry, calculated 
numerically on a grid in [9j, is therefore not accessible to this variational ansatz of 
coupled Gaussians. 

Several extensions of the Gaussian ansatz have been proposed in the literature 



10 12 , but they allow for no systematic inclusion of arbitrary angular momenta. It is 



the purpose of this work to extend the variational ansatz of coupled Gaussians in such 
a way that an instability with an arbitrary projection of the angular momentum on 
the z axis can be described for systems with axial symmetry. We present an extended 
variational ansatz including angular exponentials exp(im0), which are eigenfunctions 
of the z component of the angular momentum operator. This enables us to describe 
modes with arbitrary angular momentum within the variational framework, which will 
be demonstrated by applying the ansatz to condensates with contact interaction only, 
and with dipole-dipole interaction (DDI). 

In Sec. [2] we discuss the stability and excitations of BECs. We first present 
the method of numerically solving the BDGE, which we need for a comparison with 
the variational method. The extended variational ansatz is described and applied to 
calculate the eigenfrequencies of elementary excitations. We compare the results with 
those obtained from the full-numerical solution. In Sec. [3]we apply the variational ansatz 
to calculate the time evolution of a dipolar BEG and demonstrate that the ansatz can 
describe the angular collapse. In Sec. |4] we draw conclusions and give an outlook on 
future work. 



2. Stability and excitations 



The dime-dependence of a BEG is described in the mean-field approximation by the 
time-dependent GPE [Ij 



2M^ + y^' (p2 + AV)+iV$,„,(r,t) 



^(r,t). 



(1) 



where is the particle number, 

^Ur.t) = [ d'r'V,^,{r-r')\tP{r',t)f 



(2) 



the mean field produced by the inter-particle interaction, ipir^t) the wave function of 
the condensate normalised to unity, and M the mass of the atom species. The harmonic 
external trap is given by the trapping frequencies Up and uj^ in the p and z direction, 
and can be characterised by the trap aspect ratio A = oJzlujp- We measure all lengths 
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in units of the harmonic oscillator length aho = \/h/Mujp, all frequencies in units of uj 
and the time in units of l/up. A stationary solution of Eq. ([T]) has the time dependence 
%l){r,t) = il){r) exp(— i/xt//i), where /i is the chemical potential of the condensate. 
The interaction potential is given by [T3) 



(3) 

The quantity a is the s-wave scattering length and /im the magnetic moment of the atom 
species and is, e.g., = 6/iB (/Ub the Bohr magneton) for ^^Cr [14j. The dipoles are 
assumed to be aligned along the z axis by an external magnetic field. The interactions 
can be characterised by the dimensionless parameters Na/aho and D = N^^^ 



15 



In this article, we investigate BECs with and without LRI. For the condensate wit 
DDI we use values of D = 30 and A = 7, since for these parameters structured ground 
states appear (sjjGj. For the BEG without LRI we use a trap aspect ratio of A = -\/8, 
which was used in experiments |T6]. The scattering length can be tuned in a wide range 
via Feshbach resonances (l7j. Thus, we use this quantity as a free parameter in our 
calculations. 



2.1. Full-numerical treatment 

Besides the variational approach we also perform grid calculations to directly solve the 
GPE and BDGE, which allows us to compare both methods. Before solving the BDGE, 
one first needs to find the ground state of a BEG. In this work, we use the imaginary 
time evolution (ITE) for this task. We evolve an initial wave function on a grid in 
imaginary time r defined by r = it, and as time evolves the wave function converges 
to the ground state. The time evolution on the grid is done via the split operator 
method fls], in which the time evolution operator in imaginary time is given by 



e 



~HAr/n ^ ^-ifAr/n^-vAr/n^-ifAr/n ^ o(Ar3), (4) 



where T and V represent the kinetic and the potential terms of the mean-field 
Hamiltonian, respectively. The action of the kinetic part is trivial in Fourier space, 
whereas the action of the potential part is trivial in real space. Since we are dealing 



with a cylindrically symmetric system, we can use the Fourier-Hankel algorithm of 15 
to compute the necessary Fourier transforms for a time step and the dipolar integrals 
via the convolution theorem. 

To derive the BDGE, one starts from the usual ansatz for the perturbation of a 
stationary state ipo with chemical potential /i [Ij 

^(r, t) = [iPoir) + A (u(r) e"'"^* + v*{r) e'^*')] e-''^*/^ (5) 

where u is the frequency and A the amplitude of the perturbation (|A| ^ 1). After 
inserting this ansatz into the time-dependent GPE ([T|, while neglecting terms of 
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second order in A, and collecting terms evolving in time with exp(— iwt) and exp(ia;t), 
respectively, one obtains the BDGE 



hu)u{r) 



u[r 



(6a) 



—hjjv{r) 



v(r] 



+ J d\'Vi^t{r-r)ro{rHr')ro{r) 
+ f d\'V,^t{r-r')Mr')v{r')roir), 



(6b) 



where Hq = — + ^uj'^{p'^ + X^z'^). Due to the cylindrical symmetry, we can make 
an ansatz for the solutions 



u{r) 



u{p,z) 



vir) 



^iiruf) 



v{p,z), 



(7) 



where m is the usual quantum number of the projection of the angular momentum 
operator on the z axis. Since the Laplacian in cylindrical coordinates applied to the 
ansatz ([T]) yields 



Ae^™^M(p,2) = A 



'■p,z 



m 

7" 



u{p,z) 



(8) 



which only depends on |m|, there is a degeneracy of the eigenmodes for |m| > 0. Thus, 
we can restrict ourselves to solutions with m > 0. Another symmetry of the system is 
the reflection at the z = plane {z — )■ —z). Therefore, the solutions can be described by 
the parity vr^ = ±1 under this reflection. We denote the set of quantum numbers with 
m'^^. There are always solutions of the BDGE for w = and the trapping frequencies, 
which represent the gauge mode and the centre of mass oscillations (l]|4). 

We solve the linear, nonlocal BDGE ^ by calculating the matrix elements of the 
right-hand side in coordinate representation by means of the Fourier-Hankel algorithm 
[15j . The lowest eigenvalues of the resulting high-dimensional matrix are then calculated 
using the Arnoldi iteration |19|, implemented in the Fortran ARPACK routines [20j. 

The solutions of the BDGE ^ yield the frequencies and the shape of elementary 
oscillations of the condensates. But it may happen that the Bogoliubov spectrum of 
the ground state of a dipolar BEG contains imaginary frequencies |9], which correspond 
to dynamical instabilities. A small perturbation of the ground state then leads to an 
exponential decay. 

We calculated the ground states and Bogoliubov spectra for a BEG with an 
attractive short-range interaction (SRI) only {D = 0) and a trap aspect ration of 
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Figure 1. Bogoliubov spectra of the ground states for (a) a BEC with short-range 
interaction (A = a/8, D — Q) and (b) a dipolar BEC {\ — 1,D — 30) as a function of 
the scattering length TVa/oho (shown are the real parts of the frequencies) with even 
(solid) and odd (dashed) modes. Both spectra feature the lowest to = mode, whose 
frequency goes to zero when the scattering length reaches the critical value Ccriti ^^nd 
the centre of mass oscillations at the trapping frequencies. Unique for the dipolar BEC 
are the lowest excitations for higher angular momenta m > 0, which turn unstable and 
mark a local collapse of the condensate. 



A = -\/8j and for a dipolar BEC with A = 7 and D = 30, for which it was shown 
that the structured ground state changes its stabihty in a pitchfork bifurcation js] at 
the scattering length Na/ato ~ — 0.21£[|} The structured ground states, where the 

I In |8] a Gaussian ansatz with Nq — 5 Gaussians was used. The scattering length, where the 
bifurcation takes place, however, converges slowly with increasing number of Gaussians and may be 
different in our calculations. 
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maximum of particle density lies away from the centre, appear in isolated regions in 
parameter space with A = 7, 11, 15, . . . (6]. For our choice of parameters, A = 7 and 
D = 30, the ground state assumes a biconcave structure for the scattering lengths 
-0.546 < Na/a^o < 1.4. 

For both systems ground states only exist for scattering lengths larger than the 
critical scattering length, a > flcrit- For the condensate with SRI, the critical scattering 
length is given by Nacnt/O'ho ~ —0.457. BECs at small negative scattering lengths may 
exist, since the kinetic energy compensates the attractive interaction and stabilises the 
condensate. The critical scattering length was probed experimentally in a condensate of 



^^Rb, and significant deviations from the results of the GPE were measured 21 , which 
may be explained by taking into account higher-order nonlinear effects [22]. The critical 
scattering length of a dipolar BEC considerably depends on the trap geometry, in our 
case (A = 7, D = 30) its value is given by iVacrit/^ho ~ —0.546. A good agreement 
between experiment and the results of the GPE was found [23j. 

Fig. [T] shows the Bogoliubov spectra. Since both systems have already been 



considered in the literature (see, e. g., 24 for the BEC with SRI, and [15^j for the dipolar 
BEC), we only give a short review. The spectrum of the BEC without LRI (Fig. [T^) 
can be understood by taking the limit a — )■ 0, which yields the exact eigenvalues of the 
cylindrically symmetric harmonic oscillator 

— = 2np + \m\ + Xn^, (9) 

Up 

where rip^z ^ No are the quantum numbers for excitations along the p and z direction, 
respectively. When decreasing the scattering length towards the critical value, below 
which no stationary solutions exist anymore, the frequencies vary slightly, which means 
that in this regime the dynamics is dominated by the external trap, and the interaction 
acts as a quasi perturbation. Just slightly above the critical point the lowest mode 
with m = drops to zero, which marks the global collapse of the condensate. The 
oscillations of the centre of mass, which constantly lie at the trapping frequencies, are 
represented by the lowest odd m = mode {u/up = y/S) and the lowest even m = 1 
mode {u/up = 1). 

The Bogoliubov spectrum of the dipolar condensate (Fig. [Tjo) indicates a stable 
condensate when the scattering length is far from the critical point. However, for 
decreasing scattering length, not only the mode with m = but also frequencies for 
higher angular momenta decrease and drop to zero. At these points, those modes 
become unstable, as the eigenfrequencies are imaginary. Fig. [2] shows the real and 
imaginary parts of the lowest eigenfrequencies for the angular momenta m = 1, . . . , 5. 
The modes with m = 2 and m = 3 are the first ones which turn unstable, and thus 
mark the local collapse of the condensate. This means, that collapse occurs via local 
density-fluctuations in contrast to a global collapse, where the condensate collapse to 
the trap centre [9]. In Sec. [3] we discuss the dynamics of such a local collapse, with the 
result that an extended variational ansatz is necessary to describe the collapse. This is 



A variational approach to Bogoliubov excitations and dynamics of dipolar BECs 



7 



3 
3 



3 
3 



1.4 
1.2 
1 

0.8 
0.6 
0.4 
0.2 


0.2 
0.4 
0.6 
0.8 

1 

1.2 
1.4 



-0.6 



-0.4 



-0.2 



m=l, real 
m=2, real 
m=3, real 
m=4, real 
m=5, real 





Na/ a, 



m=l, imag. 
m=2, imag. 
m=3, imag. 
m=4, imag. 
m=5, imag. 

0.2 



0.4 



0.6 



ho 



Figure 2. The lowest modes of the Bogohubov spectrum of the dipolar BEC from 
Fig. [TJd. Shown are the real and also imaginary parts of the frequencies. In our 
resolution of the scattering length, the modes with m = 2 and m = 3 turn unstable 
simultaneously, but the latter has a higher imaginary frequency below the critical 
scattering length, which means that this mode dominates the collapse in the vicinity 
of the stationary solution. 



a unique feature of dipolar BECs and is related to the biconcave structure of the ground 
state occurring for specific trapping geometries (8|[9]. These instabilities are referred 
to as rotons, or, more specifically, as discrete "angular rotons", since these excitations 
assume a shape with azimuthal nodal surfaces [6]. 

Worth mentioning also is the avoided crossing in the lowest modes for the angular 
momenta m = 0, 2, 3, whereas for m = 1 there is a regular crossing since u = Up is an 
exact solution for every scattering length, which does not allow for an avoided crossing. 



2.2. Variational approach 

The grid calculations are very accurate, if grid size and resolution are chosen carefully, 
but computationally expensive. The variational method provides an alternative in that 
the computation time reduces significantly compared to grid calculations. With this 
motivation, Rau et al (2||3|[8] used the variational ansatz 

Ng 
k=l 

to calculate stationary solutions and excitations of dipolar condensates. Here, Nq is 
the number of Gaussians, A^, Ay and A^ describe the widths and 'y^ the amplitude 
of each Gaussian. The quantities {A^, Ay, A'^,'y'') are the variational parameters of 
this ansatz. In a recent article [4j we showed that the variational ansatz with coupled 
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Gaussians can only describe modes with an angular momentum of / < 2, \m\ < 2. 
For that reason we proposed an extended variational ansatz with coupled Gaussians 
and spherical harmonics, which can describe modes with arbitrary angular momenta in 
spherically symmetric systems. We found a good agreement of the lowest eigenmodes 
with those obtained by grid calculations. 

For the systems with axial symmetry considered here, this ansatz is not applicable, 
since the angular momentum is no longer a good quantum number and therefore the 
angular part of the excitations is not described by spherical harmonics. However, the 
z component of the angular momentum, represented by the quantum number m and 
the eigenfunctions exp(im0), is still conserved. Using this knowledge we construct an 
extended variational ansatz, which is specific to cylindrically symmetric systems, and 
which is given by 

^ = Y^il + ^Y^ 4_ppl™U^'e^"^M e-^>'-^5^'-P'^-^'. (11) 

k=l \ m^Op=0,l / 

The complex variational parameters and describe the width in the p and z 
direction, is the displacement along the z axis, and describes the amplitude of 
each Gaussian. With these parameters alone, the wave function does not depend on 
the angular coordinate 0, thus the parameters {A^, A'^,p^,'y'') are responsible only for 
any m = excitation. To account for higher angular momenta and arbitrary parity 
TTz, we include the sums over m and p in the brackets (vr^ = +1 belongs to p = 0, and 
TTz = —1 to p = 1). The amplitude of each of these excitations is given by the variational 
parameters d^p. 

To derive the equations of motion for the variational parameters, we follow the 



procedure in |2||4j and apply the Dirac-Frankel-McLachlan TDVP 25 , 26 . It states that 
the norm of the difference between the left- and the right-hand side of the Schrodinger 
or GPE 

i = \Mt)-Hm\\' (12) 

must be minimised for a given variational ansatz ip = ip{z) with respect to the variational 
parameters z. The variation of while ip is fixed, and the identification of x with ip 
leads to the necessary condition (27| 



Kz = -ih, (13) 
where the matrix K and the vector h are defined by 

H^p). (14) 



dipX ^ / dip 

dzj / ' ' \dzi 



The remaining task is to evaluate all integrals appearing in Eq. (14). Technical details 
are presented in Appendix A Since the ground state in the GPE corresponds to a fixed 
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point of the equations of motion (13), we can find the ground state of a system by a 



nonhnear root search, for which we require 



i/i for Zi = 7* 
else, 



(15) 



where fi is the chemical potential. 

Stability and excitations in the framework of the TDVP are given by the linearised 
equations of motion, for which one first needs to rewrite the complex equations of 



motion (13) into real ones by defining the new vector of variational parameters z as 



z = (Rez, Imz). The time-dependence of a small perturbation Sz is then given by ||2j 



62i = J2 



dzj 
dz, 



(16) 



z=zo 



with the Jacobian J evaluated at the fixed point Zq. These linearised equations of motion 
are as usual solved by a harmonic time-dependence of the perturbation 



6z{t) 

which leads to the eigenvalue problem 



'lUlt r ~ 



5i(0), 



iuSziO) =JSz{0). 



[17) 



(18) 



A real u belongs to an oscillatory motion, as for real eigenvalues u of the BDGE. 
The symmetry of the system, which leads to the quantum numbers m and vr^ in the 
Bogoliubov spectrum, gives rise to a block structure of the Jacobian which is further 



explained in [Appendix B 



The eigenvalues of the Jacobian and the BDGE can be directly compared. But 
is there also a relationship between the eigenvectors of the Jacobian 6z and the 
eigenf unctions of the BDGE u and For a given eigenvector 6z, a wave function can 
be constructed by taking the difference of the excitation with variational parameters 
Zq + Sz and the wave function of the stationary solution with the parameters Zq, 



Sipizo, Sz) = ij{zo + 6z) - tp^zo) 



(19) 



Since we are dealing with an excitation in the linear regime, one can Taylor expand the 
first summand in Eq. (19) to first order in Sz. This yields the function |2] 



Sip = Sz 



dip 
dz 



(20) 



=Zo 



which is the scalar product of the eigenvector and the gradient of the wave function 
with respect to the variational parameters evaluated at the stationary solution. Thus, 
every eigenvector Sz belongs uniquely to a function Sip. 
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As the Jacobian J is real and a stable mode is represented by an imaginary 
eigenvalue icu, for every eigenvalue iuj with w > there is also an eigenvalue — iw, 
therefore we can construct corresponding functions 5ip+ and 5ip- which evolve in time 
as exp(ici;t) and exp(— iwt), respectively. They can be combined to the perturbation 

5^ = 5^_e-'^' + 5ij+e''^\ (21) 

Comparing this with the Bogoliubov ansatz from Eq. ([s]) leads to the relationship 

u^5i)-, v*^5il}+, (22) 

which links the eigenvectors of the Jacobian with the eigenf unctions of the BDGE. 



2. 3. Results for the EEC with scattering interaction 



We now apply the variational ansatz (11) to a condensate with an attractive SRI. 



Fig. [3] shows the results, where A^g = 5 coupled Gaussians have been used (solid lines), 
compared to the full-numerical calculations (dotted lines). For the sake of clarity we 
distinguish between even modes in Fig. |3^ and odd modes in Fig. |3)d. One recognises 
that the lowest mode of each angular momentum and parity can perfectly be reproduced 
by the variational ansatz. Interestingly, as with our previous ansatz for spherically 
symmetric systems |4|, the frequencies of the centre of mass oscillations (lowest odd 
m = at u; = -\/8wp, lowest even m = 1 at a; = UJp) agree, within numerical accuracy, 
with an arbitrary number of coupled Gaussians, even with A^g = 1- Foi' the second 
lowest modes we also find a very good agreement as long as the scattering length is not 
close to the critical value. For the case m = 0, vr^ = 1 even more eigenvalues of the 
Jacobian are close to the numerical results. 

To demonstrate the power of our extended ansatz, we calculated eigenmodes for 
angular momenta up to m = 6 (see Fig. |4]) for a fixed scattering length Na/aho = 
—0.34. Even for these angular momenta the lowest modes agree well. For the second 
excitations and lower angular momenta we find also a good agreement, however, the 
deviations of the variational and numerical results become larger for increasing angular 
momenta. Nevertheless, there is an obvious relationship between the eigenfrequencies 
of the Bogoliubov spectrum and the eigenvalues of the Jacobian, where our extended 



variational ansatz (11) - specific to cylindrically symmetric systems - has been used. 
Thus, the variational ansatz is a valid alternative to simulations on grids if one is 
interested only in the lowest modes. 



2.4- Results for the dipolar EEC 

We now want to include the DDI in the calculations and especially verify if the extended 
variational ansatz can describe the stability change for higher angular momenta m > 0. 
We first investigate the eigenmodes in the stable regime with a scattering length 
A^a/oho > 5. Fig. Is] shows the comparison of both spectra for the angular momenta 
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Figure 3. Comparison of the Bogoliubov spectrum for a BEG with SRI from Fig. [iji 
with the eigenvalues of the Jacobian obtained from the variational ansatz for even (a) 
and odd (b) parity. We find a very good agreement for the lowest two modes of each 
angular momentum and parity, only close to tlie critical scattering length deviations 
of the second lowest modes for to > can be seen, which become larger for increasing 
angular momentum. 



m = 0, . . . , 3. With 6 coupled Gaussians, the lowest modes of each angular momentum 
and parity (solid lines) are converged to their corresponding numerical frequencies 
(dotted lines). As in the case of the BEG without LRI, the centre of mass oscillations 
agree within numerical accuracy. For the modes with even parity, there is a good 
agreement even for the second lowest modes, where the differences become larger for 
rising scattering length and angular momentum. In the case of odd parity, the deviations 
for the second lowest modes become quite large for m > 2, when the scattering length 
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Figure 4. Comparison of both spectra as in Fig. [3] but here for a fixed scattering 
length of Na/a\io — —0.34 and angular momenta up to m = 6. The variational 
ansatz also works for these modes and can describe the lowest modes of the Bogoliubov 
spectrum, therefore we are able to establish a connection between both methods. 



is large. 

We have also investigated eigenmodes with an angular momentum up to m = 6 in 
the stable regime. Fig. [6] shows the results for a fixed scattering length Na/aho = 6. 
The variational ansatz can reproduce the lowest modes even for these angular momenta 
and, since the scattering length is not too large, also the second lowest modes. For the 
case m = 0, TTz = I, we find the lowest three modes converged, and a good agreement 
of the fourth lowest mode. These investigations show that the variational ansatz (11) 
can also reproduce the full-numerical Bogoliubov spectrum, if only the lowest modes are 
compared. 

We now investigate the interesting case of a small negative scattering length, where 
the Bogoliubov spectrum shows unstable roton modes for different angular momenta 
(see Fig. [2|. We focus here on the case m = 3, which cannot be described by the ansatz 



of coupled Gaussians Eq. (10), hence the extended variational ansatz is necessary. Fig. [7] 
shows the lowest modes with m = 3 obtained using different numbers of Gaussians (solid 
and dashed lines). For comparison the full-numerical results are also plotted (dotted 
line) . We find that with an increasing number of coupled Gaussians the accuracy of the 
variational calculations increases gradually. 

Other angular momenta show a similar behaviour, the convergence is in general 
better for lower angular momenta. Thus we can conclude, that the variational 
ansatz (11) can describe, though not quantitatively perfect, the dynamical instabilities 
of the Bogoliubov spectrum for angular momenta m > 0. 

It is known that at the critical scattering length below which no stationary solution 
exists anymore an unstable collectively excited state emerges in a tangent bifurcation J3] 
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Figure 5. Comparison of the Bogoliubov spectrum for a dipolar BEG from Fig. [TJd 
with the eigenvalues of the Jacobian obtained from the variational ansatz. Shown are 
the modes with even (a) and odd (b) parity in the stable regime witli a scattering 
length A^a/flho > 5. Whereas the two lowest modes agree well for even parity, for 
odd parity the second lowest modes differ quite clearly, especially for higher angular 
momenta m > 2. 



in addition to the ground state (see the inset in Fig. [8]). For condensates without DDI, 
this excited state is unstable with respect to an m = excitation [s]. The variational 
approach offers the possibihty to investigate this state of a dipolar BEC, which is not 
accessible to the full-numerical ITE method, and its collapse mechanism. Fig. |8] shows 
the imaginary frequencies of the excited state obtained with Nq = 5 Gaussians. As both 
states merge at the critical scattering length the same holds for the eigenfrequencies, 
and the excited state shares the instabilities of the ground state at the critical scattering 
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Figure 6. Comparison of both spectra as in Fig. [5] but here for a fixed scattering 
length of Na/a\io = 6. For the lowest mode of each angular momentum and parity, 
there is a very good agreement. For the second lowest modes we find differences, which 
get larger for increasing angular momenta. The correspondence is slightly better for 
even than for odd modes. As in the case without LRI there is a link between the 
full-numerical Bogoliubov spectrum and the eigenvalues of the Jacobian resulting from 
the variational ansatz. 
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Figure 7. Lowest mode with to = 3 which turns unstable when the scattering length 
is decreased. Compared are the full-numerical frequencies (dotted line) with the results 
from the variational ansatz with different number of Gaussians (solid and dashed lines). 
The accuracy of the variational calculations increases gradually with the number of 
coupled Gaussians. 
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Figure 8. Imaginary frequencies of the excited state with Nq = 5 Gaussians. At the 
critical point the excited state shares the same instabilities as the ground state with 
additionally an rn = instability, which is typical for the collectively excited state. The 
inset shows the tangent bifurcation in the mean field energy i?i„t/7V£^ho {N being the 
particle number) for the ground (solid) and excited (dashed) state, where E'ho = 
is the energy quantum of the harmonic oscillator. 



length, except the additional m = instability. In contrast to the case with SRI only or 
with monopolar LRI Q the excited state of the dipolar BEC for the trapping frequencies 
is unlikely to collapse with m = symmetry due to the additional instabilities for m > 0, 
which have a higher imaginary frequency and thus a larger collapse rate. 

The m = 5 excitation in Fig. |8] changes its stability with increasing scattering 
length, in a similar way as the unstable modes of the ground state (see Fig. [?]). Note, 
however, that this does not mean a global change of the stability due to the other 
unstable modes. The excited state stays unstable for all scattering lengths considered 
in our calculations. 



2. 5. Shape of Bogoliubov excitations 

After linking the eigenvalues of the BDGE and the Jacobian of the time-dependent 
variational approach, we will check if there is also a direct correspondence between the 



eigenfunctions and eigenvectors of these two approaches. In Sec. |2.2| we discussed how 
to construct the functions dip- and dip^ from the eigenvectors and how they are related 
to the functions u and v. For a better comparison of both methods, we use the usual 
normalisation of the Bogoliubov functions J d^r \\u'^ — \v\^] = 1 jl], and assume the 
same for the variational method, J d^r[\6il!-f — \6ip^f] = 1. We point out that in 
calculating 6il>± there is no ambiguity, except for a global phase, which we choose in 
such a way that the functions are real. Fig. [9] shows the comparison of both methods 
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Figure 9. Radial dependence for the functions u, v and Sip^, 5ip+ (red lines), 
respectively for different values of z [0 (blue crosses), 0.42aho (yellow squares), 0.84aho 
(green circles)], where A^g = 8 Gaussians have been used (see text for normalisation). 
Shown are the first two excitations for the angular momenta m = 0, 3 and a fixed 
scattering length of Na/uho = 6. All functions are chosen to be real. There is a very 
good agreement, only for the second m = 3 excitation there are small deviations. 



for some modes. The first two functions u with m = have one and two radial nodes, 
respectively, because the lowest mode with m = represents the gauge mode, which is 
no physical excitation. For other angular momenta, the n-th excitation in p direction 
has n — 1 radial nodes. The asymptotic behaviour of all functions in Fig. |9] for p — )■ 
is p'™"'. The nodal structure for the functions v is different. All f for m = have one 
node, and for the other angular momenta zero nodes. For higher excitation numbers 
and angular momenta, the amplitude of v becomes smaller and smaller, indicating that 
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in these limits the couphng between u and v in the BDGE can be neglected, which is 
in accordance with previous calculations in the literature [l, 15,28). 



For the second excitation with m = 3 we find small deviations, which is to 
be expected, since there are also deviations in the eigenfrequencies (cf. Fig. [6]). By 
comparing more eigenmodes we can conclude: The better the agreement of the 
eigenfunctions, the better is also the accordance of the frequencies. For higher 
modes, where there is no quantitative agreement between the eigenfrequencies, only 
the asymptotic behaviour and the nodal structure can be described by the variational 
ansatz, even though there is at least a qualitative agreement in those aspects. 

The comparison of the eigenfunctions of the BDGE with those obtained from the 
eigenvectors of the Jacobian reveals a very good agreement. This is remarkable, since 
both functions originate from quite different approaches of stability. The extended 



variational ansatz (11 ) cannot only quantitatively describe the frequencies of the lowest 
modes, but also the shape of the oscillations, even for angular momenta, which were not 
accessible by the variational approach without the extension. 

3. Dynamics 



The extended variational ansatz (11 ) was proposed to find a direct relationship between 
the eigenvalues of the Jacobian and the eigenfrequencies of the BDGE. We have derived 
the equations of motion using the TDVP to calculate the Jacobian. The equations of 



motion (13) are not only valid in the vicinity of a fixed point, but for all variational 
parameters for which the integrals in the equations of motion converge. Wilson et 
al |9) simulated the dynamical evolution of a collapse for a dipolar condensate with 
biconcave structure by reducing the scattering length below the critical value, which 
yields a symmetry breaking collapse {angular collapse) with m = 3 symmetry. Rau 
et al [3] simulated this scenario within the variational framework with an ansatz of 



coupled Gaussians. However, as mentioned in Sec. 2^, this ansatz cannot describe 
wave functions that exhibit the m = 3 symmetry. Motivated by these papers we will 
now integrate the equations of motion numerically with appropriate initial conditions to 
verify, whether or not the extended variational ansatz is also applicable in calculating 
the dynamics of dipolar condensates and especially to describe the angular collapse 



for arbitrary angular momenta, which is caused by the roton instability (cf. Sees. 2.1 



and 2.4). To do so we only address the collapse of ground states with biconcave structure. 

As a first example, we consider the system with Nq = 4 Gaussians and calculate 
the dynamics of an induced collapse by lowering the scattering length below the critical 
value. As an initial state for t = 0, we take the stable ground state for a scattering 
length of Na/aho = 6 and add a small perturbation for a specific angular momentum 
by changing the variational parameter d^^^Q to a small value larger than zero, e.g. 



c/^p^o = 2 X 10 ^, which is just a small perturbation of the ground state. We let 



•,,p= 

m,p=0 

the condensate evolve until tUp = 0.72, after which the scattering length is linearly 
ramped down until tUp = 10.8 to a value of Na/a^^o = —0.6, which is below the critical 
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(a) m=3, trap=20.88 (b) m=4, trap=25.99 




Figure 10. Modulus squared of the wave function of the collapsing dipolar 
condensate in a cut through the [z — 0)-plane. Bright (yellow) marks regions of 
high particle density, while dark (black) stands for low densities. For the calculations 
Nq = 4 Gaussians were used. 



scattering length. The condensate then evolves at a fixed scattering length, until it 
collapses to a point, from where on no further integration is possible. We calculated 
this scenario for an initial m = 3 and m = 4 excitation, which cannot be described 
by the variational ansatz without the extension and thus is an advance as compared to 
previous calculations of a collapsing condensate within the variational framework [3]. 



Fig. 10 shows the resulting density of the collapsing condensate for both cases. For 
an initial m = 3 excitation, the collapsing cloud exhibits the clear shape with an m = 3 
symmetry, the angular dependency is proportional to sin(30). This result shows that the 
extended variational ansatz is not only able to describe the instability in the vicinity of 



the ground state, but also the nonlinear dynamics given by the equations of motion (13). 
If the initial state is perturbed by an m = 4 excitation, the resulting density distribution 
shows a superposition of several angular momenta. There are altogether four density 
peaks, of which two are significantly larger than the other two. The appearance of four 
peaks corresponds to an m = 4 symmetry, but the distribution of the particle density 
clearly belongs to m = 2. Of the two large peaks, one is slightly higher than the other 
one which indicates an m = 1 symmetry. The calculation shows that the different 
angular momenta do not evolve freely, but instead are coupled and infiuence each other. 
This is due to the nonlinearity of the GPE ([T]) which carries over to the equations of 
motion for the variational parameters. 

There are many excitations one can add to the ground state as a perturbation. 
Above, we added specific eigenmodes to investigate the behaviour of the extended 
variational ansatz when considering the collapse dynamics. We now want to utilise a 
more physical motivated initial state. In Ref. |9] the authors added modes with weights 
determined by the Bose-Einstein distribution, which leads to 



where ifjQ is the ground state. The weight of each excitation, where m is the angular 
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Figure 11. Particle densities of the collapsing cloud at tujp — 12.384 with high (low) 
density represented by bright (dark) regions. All calculations were done with Nq = 5 
Gaussians. The different random phases used to construct the initial state (see text) 
are the only differences in the four calculations. There are three density peaks each, 
but the overall rotation and the density distribution along the peaks are different. 



momentum and n summarises the other quantum numbers, is given by 



(24) 



Here, 



is the energy of each excitation, k-Q the Bohzmann constant, and T the 



temperature. The quantities an,m in Eq. (23) determine the initial phase of each 



excitation and are random numbers between and 1. The overall factor l/y/N, with 
the particle number, is necessary, since the wave function ip is normalised to unity 
instead of A^. 

We now can make use of the relationship between the eigenvectors of the Jacobian 



and the eigenfunctions of the BDGE (Sec. 2.5) to prepare an initial state similar to that 



given by Eq. (23) within the variational framework, where instead of the functions, the 



eigenvectors of the variational parameters have to be added to the ground state vector. 
For the following calculations we simulate a condensate consisting of = 10000 ^^Cr 
atoms at a temperature of T = 100 nK. We start with an initial wave function according 
to Eq. (23) and a scattering length of Na/a^io = 6, which is ramped down to a value of 
A^a/ttho = —3 over a time of AtWp = 14.68. 

We performed several calculations with this scenario, for which the only differences 
are the initial phases an,m- As the scattering length decreases the particle density 
establishes its biconcave structure, which is a typical feature of dipolar BECs (6]-[9). If 
the scattering length further decreases, the critical value is reached and the atomic cloud 
starts to collapse to the region with the highest density, which is on a ring around the 
centre. Due to the initial excitations with angular momentum m > 0, which are now 



unstable, the symmetry is broken and several density peaks emerge. Fig. 11 shows 
the density distribution of four different calculations shortly before collapse. Each 
condensate has three peaks located almost equidistantly on a ring. But the distribution 



along the ring differs: In Figs. 11 (a)-(b) there are two large peaks and a smaller one, 



while in Figs. 11 (c)-(d) it is the other way round. Additionally the initial phases 
determine the overall rotation of the collapsed condensate in the {z = 0)-plane. We can 
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conclude that both, the m = 2 and m = 3 mode, are responsible for the shape of the 
collapsed cloud, which corresponds to the numerical solution of the BDGE (see Fig. [2]), 
in which both modes turn unstable almost simultaneously. 

A question, which naturally arises, is how the number of Gaussians influences the 
dynamics calculation. We performed similar calculations with Nq = 4 and Nq = 6 
Gaussians. With 4 Gaussians we find a different behaviour in that a collapse with 
m = 2 symmetry is favoured. However, with 6 Gaussians almost the same results as in 
the case with 5 Gaussians are obtained. The only difference is the time taken before the 
condensate collapses, which is shorter with 6 Gaussians. This can be explained by the 
fact that the critical scattering length is shifted to a higher value and thus the critical 
point is reached earlier during the simulation. The eigenfrequencies with our Gaussian 
ansatz do not fully agree with the full-numerical results (cf. Sec. 2.4), however, we expect 
the variational results to not being completely converged with 5 or 6 Gaussians. 



4. Conclusion and outlook 



We presented an extended variational ansatz, which is specific to cylindrically symmetric 
systems and which can describe arbitrary angular momenta, and derived the equations 
of motion for the variational parameters using the TDVP. The eigenfrequencies of two 
different systems (without and with DDI) were calculated using the variational ansatz 
and compared with the direct solutions of the BDGE. We found a good agreement of 
the lowest-lying modes and could confirm the direct relationship of both methods. This 
extended our investigations of spherically symmetric systems |4] and we could show that 
this agreement is of generic type, independent of the symmetry and the interactions. 

The roton-instability of a pancake-shaped dipolar BEG was investigated using 
the extended ansatz. We could find instabilities for arbitrary values of the angular 
momenta, which is a great progress in comparison with the variational ansatz with 
coupled Gaussians without the extension. With more and more Gaussians we obtained 
gradually more accurate results within the variational framework. Thus, our variational 
method is a valid approximation to grid calculations. 

We further constructed functions of the eigenvectors of the Jacobian and found 
that these functions agree very well with the solutions of the BDGE, as long as 
the corresponding eigenvalues agree. There is not only a direct link between the 
eigenfrequencies, but also between the eigenvectors and the eigenf unctions. 

After considering the dynamics in the linear regime we directly integrated the 
nonlinear equations of motion and showed that our variational ansatz, assuming that 
the initial conditions and the scenario of the time-dependent scattering length are 
appropriately chosen, is capable of describing the angular collapse with arbitrary angular 
momenta. We used the relationship between the eigenvectors of the Jacobian and 
the eigenfunctions of the BDGE to construct an initial state which corresponds to a 
condensate in thermal equilibrium within the mean field approximation, and simulated 
a realistic experiment of the angular collapse. We found that the random initial phase 
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distribution strongly influences the shape of the collapsed atomic cloud, which makes 
the resulting shape of an experiment irreproducible. 

Altogether it was shown that the extended variational ansatz is a significant progress 
compared to coupled Gaussians, and as yet unaccessible results could be obtained. 
Our ansatz is a valid alternative to calculate eigenfrequencies and shapes of elementary 
excitations, if the lowest eigenmodes of stable condensates are considered. The roton- 
instabilities can be obtained, but we could not find a good quantitative agreement, 
leaving our ansatz an approximation. 

In future applications of the extended variational ansatz it would be possible to 
include the particle loss in the dynamics calculation, as the losses are necessary to 
accurately simulate a collapse scenario of a dipolar BEG (91129). Furthermore, our ansatz 



is applicable to dipolar BEGs in a one-dimensional optical lattice 30 32 to calculate the 



ground state, stability, excitations and the dynamics. The direct relationship between 
the eigenvectors of the Jacobian and the eigenfunctions of the BDGE allows several 
applications of the variational approach, e. g. to describe a BEG at finite temperature, 
where a self-consistent solution of the GPE and BDGE is necessary l33). 
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Appendix A. Calculation of the integrals for the variational ansatz 



This appendix is dedicated to the integrals necessary for the matrix K and the vector 
which are defined by 



dip 
dzi 
dip 
dzi 



dijj 
Hijj 



Hip{r) 



(A.la) 
(A.lb) 



to set up the equations of motion (13) of the variational ansatz (11) via the time- 



dependent variational principle. We write the variational ansatz in the form 



Ng 



m,pr ^ c c , 



(A.2) 



k=l m p=0,l 



where d^Q = 1 and = are no variational parameters, but constants. We further 
introduce the abbreviations 



{Zi\Zj) = 



dijj 
dzi 



dip 

dZj 



(z^lXi/j) ^ 



dip 
dzi 



Xip 



(A.3) 



for the matrix elements, where X is an operator of the mean field Hamiltonian 
H = T + I4xt + $s + '^'d with the kinetic part, the external potential and the scattering 
and dipolar mean field potentials. 
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We do not give a detailed way on how to calculate all necessary integrals, nor do we 
present all results. We rather want to sketch the calculations and give results of some 
examples. For all details we refer to the master thesis of one of the authors |34| . 

Appendix A.l. Integrals of the K-matrix 

Before giving the results, we define the abbreviations 

Aj'^Aj + (4r, Af^A^ + (Air, + ^'^^ ^ + (V)* (A.4) 

and the basic integrals 

oo oo 







The p integral can be identified with the gamma function by a simple substitution. 
The z integral is a Gaussian integral for n = 0, higher orders can be obtained by the 
recursion formula 

d_ 
dpi 

The integrals of the matrix K then yield 



(Cp2l^mi,pi> = 2Ti5m,mj''\\mi\ + 1,^1+^2), (A.7a) 
(CpJ4> = -2TiY.l'\\m^\+2,p,+p^), (A.7b) 



pi 



dl.J^l) = -27r5^/'='(|m2| + l,Pi+P2 + 2), (A.Tc) 

'A.,vM) = -27r5^/'^^(|m2| + +P2 + 1), (A.7d) 

pi 

;Cpj7'> = -2^Y.^^\\m,\ + (A.7e) 



pi 



The other integrals can similarly be expressed by Eq. (A.5). 



Appendix A. 2. Integrals of the kinetic term 

For the integrals of the kinetic term the action of the Laplacian on the variational ansatz 
is necessary. This is most easily evaluated in cylindrical coordinates. The integrals can 
then be expressed by Eq. ( |A.5[ ). For the d parameters, e. g., one obtains 

<Cp. I^P = ^ E E { [4(|m2| + + 2{2p, + l)Al - [p^f] 

k pi 

X I^\\m2\ + +P2) + [2pipf\ /''(|m2| + 2,p2) - [4^^^^ /''(hsl + 2,pi +p2 + 1) 
- [4(A^r]/''(KI+2,pi+p2 + 2)- [4(Aj)2]/^''(|m2|+3,pi+p2)|. (A.8) 
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Appendix A. 3. Integrals of the trapping potential 

The trapping potential is of the form \4xt/ = ^^I^P^ + A^z^). Thus, the corresponding 
integrals are easy to evaluate. As an example we obtain 



k pi 

X (^J^'(|m2| + 2,pi+p2) + A2/^-'(|m2| + l,Pi+P2 + 2)). (A.9) 
Appendix A. 4- Integrals of the scattering interaction 

The integrals for the interaction terms require some more attention. First an expression 
for the modulus squared wave function {ipf is needed. This enters the integrals for the 
interaction potentials. The scattering interaction is easier to calculate, since this kind 
of interaction is local, in contrast to the DDI. 

We further need new abbreviations, which depend on four indices and which are 
defined by 

^ = A;^' + Aj', = A^^' + Af, pT^=pi'+pf, 7*^''=' =7^^' + 7'='. (A.IO) 



The new basic integral, P^''\m,n), can be obtained from the former one, Eq. (A.5), by 
substituting all indices kl with ijkl. The scattering integral then reads 

(^m2,P2 l^sV') /^p = 8^^- ^ ^ ^ ('^m4,P4)*'^m3,P3^mi,pi'^mi+m3,m2+m4 

Who . . , 

i,j,k mi,r?i3,m4 pi,p3,p4 

XP'^'r ^-^ ,Pl+P2+P3+Pij ■ (A.ll) 

Appendix A.5. Integrals of the dipolar interaction 

In order to calculate the integrals of the dipolar interaction, we first need the expression 
for the mean field potential $d- This potential can be rewritten using Fourier transforms 
and the convolution theorem, which yields (35) 



Mr) = J^-'[VM-Hk)], (A.12) 
where the functions with tilde designate the Fourier transforms of the dipole potential 



and the particle density. The transform of the dipole potential is well-known to be 35 



47r / A;2 \ 

V,{k)/nw, = - (^3^ - ij D/N. (A.13) 

Our next task is to calculate n{k). Writing the square of the absolute value of the wave 
function with the given variational ansatz, one notices that the angular dependency is 



of the form exp(im0) with some m G Z. In 15 it was noted that the three-dimensional 



Fourier transform of such a function can be expressed as a one-dimensional Hankel 
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transform for the p direction and a remaining one-dimensional Fourier transform in the 
z direction. The Hankel transform brings in the Bessel function, and the integral can 
be worked out by Taylor expanding the Bessel function. The Fourier transform is a 
Gaussians integral with linear terms and variable powers of z, for which an analytical 



expression can be found in established tables of integrals 36, Eq. 3.462.2]. 

In the original integral, ((i|<l>dV')) the integral over r is similar to the Fourier 
transform of \'(p\^ and can easily be expressed in the same way, which yields the 
intermediate result 

i,j,k mi,m3,m4 pi,p3,p4 
•mi— r7i2— m3+m4 

y ( Akl\-\ml-m2\-^lml.m2-'^( A^j)-\"^3-mi\-fi„i^^„i^-l 

2|mi-m2| + |m3-m4|+2'^ P'' l^p J 



X ^ ^ j ^^rn^<'^'2 \ j ^^'rn3,m4 \ ^^^fc«^Qj^^^ij^/3^_]^^Mmi,m2+A""3.'"4 
X (- |mi - m2| - ^mi,m2)M™i,™2-" l"^3 ~ "^^1 - l^rm^rm) t,^^,^^-^ X Jfc, (A.14) 



where we defined 

a ^fsimmr u =1° for sign ^ sign m^, 

I mm(|mi| , \m2\) otherwise, 

and used the Pochhammer symbol (a)n, which is defined by 

(a)„ = a(a - l)(a - 2) • ■ ■ (a - n + 1) = (A.16) 

1 (a — n + 1) 

The remaining three dimensional integral 1^ reads 



with the z integrals 



oo 

rj{p)= j dzzPe-^''''-^^'+'^'^\ (A.18a) 

— oo 

oo 

/!',(p)= j dzzPe-'^^'^'-^P''-'''^^', (A.18b) 



which are evaluated in Ref. 36, Eq. 3.462.2]. The Fourier transform of the dipole 



potential consists of two parts, Vd = {Airk'^/k'^ — 47T/3)hu)pD/N = The integral 
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of Ik containing is equivalent to the integrals of the scattering interactionjH and 
requires no further attention. The part of with V^, which we denote by can be 
rewritten to 

Il/hujp = 47r2(D/iV)5„,+„3,„,+^,n!(af 0" 

oo oo 

X y"dtt-("+i) j dhkle-^'-^^''^''''i'J,{p,+p2)lHp3+P4). (A.19) 

1 — oo 

The kz integral is of Gaussian type, and similar to J^'^ and Pj . The remaining t integral 
can be expressed by the hypergeometric Function 37 , assuming all parameters 
are zero. In the general case, the t integral can be efficiently evaluated using a Taylor 
series. 



Appendix B. Structure of the Jacobian 



The solutions of the EDGE Q can be classified by the quantum numbers m and vr^, 
which express the symmetries of the system. Applying the TDVP and linearizing the 
equations of motion yields an alternative description of linear oscillations and stability, 
in which the Jacobian J plays a central role. The Jacobian, similar to the EDGE, also 
express the symmetries of the system and the quantum numbers, in that it assumes a 
block structure, which can in general be written as 



{ J±m} 



(B.l) 



•7 



where is the submatrix belonging to the angular momenta ±m and the parity 

TTz- This block structure allows one to set up and diagonalise the Jacobian for each 
m and vr^ individually, analogously to solving the EDGE. The only difference is that 
in the Jacobian the angular momenta m and — m are coupled, which is due to the 
coupling of those angular momenta in the integrals necessary for the equations of motion 



(cf. Appendix A) 



A simple example, which shows that the coupling may not be neglected, is the 
submatrix of the m = ±1, vr^ = 1 Jacobian obtained with A^^g = 1 Gaussian, for a 
dipolar EEC with A = 7 and D = 30 



/ 

2.5936 


y-2.3931 



-2.5936 


-2.3931 






-2.3931 


2.5936 



-2.393l\ 


-2.5936 
/ 



pi 



(B.2) 



§ The Fourier transform of 14 = 47riVa/aijo(5(r) is AirNaj a\io, and thus proportional to V^. This means 
that the whole integral of this part is equivalent to the integrals of the scattering interaction. 
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where the variational parameters are ordered as 2 = (Re (i^ q, Imd^ g, Re rf'^^ q, Imd^ q) 
for a scattering length of Na/a-^o = 6. The eigenvalues of the matrix (B.2) are 
(+iWp, +iu}p, —iUp, —iujp). The eigenvectors of the eigenvalue +10;^ can then be combined 
to an eigenvector belonging to m = 1 and another eigenvector for m = — 1, the same 
holds for the eigenvalue —iUp. In general, for arbitrary number of Gaussians and angular 
momenta, by taking appropriate linear combinations of the eigenvectors, eigenvectors 
belonging to a specific angular momentum +m or —m can be obtained. 
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